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Abstract 

Convergence properties of a multigrid algorithm, developed to calculate compressible viscous flows, are 
analysed by a vector sequence eigenvalue estimate. The full three-dimensional Reynolds- averaged Navier- 
Stokes equations are integrated by an implicit multigrid scheme while a k — turbulence model is solved, 
uncoupled from the flow equations. Estimates of the eigenvalue structure for both single and multigrid 
calculations are compared in an attempt to analyse the process as well as the results of the multigrid 
technique. The flow through an annular turbine is used to illustrate the scheme’s ability to calculate complex 
three-dimensional flows. 


Introduction 

In many ways the use of computational fluid dynamics has been legitimised by the successful calculation 
of two-dimensional inviscid transonic flows. However, for computational fluid dynamics to have a continued 
impact on the aerodynamic design process, efficient three-dimensional methods must be developed. Because 
the greatest aerodynamic gains are now being realised by truly three-dimensional designs, the need to 
compute three-dimensional compressible viscous flows is apparent. 

The present work deals with the analysis of a multignd scheme to solve the Reynolds- averaged Navier- 
Stokes equations. The flow equations are solved by the diagonally inverted LU implicit multigrid scheme , 1 
while the Reynolds stess tensor is modelled by a Boussinesq eddy-viscosity formulation which calculates 
turbulent viscosities from a high Reynolds number k — e model. 

The multigrid method has been, by far, the most successful means of accelerating inviscid calculations 
to steady state, A number of researchers have recently extended this technique to viscous flow calculations. 
Although it is true that relatively greater benefits are produced with the Euler equations, a significant 
amount of acceleration is also produced with the Navier-Stokes equations. In an attempt to understand the 
process by which the multigrid method produces this acceleration, eigenvalue estimates are compared for 
both single and multigrid calculations of inviscid and viscous flows. 


Analysis 

The three-dimensional Reynolds-averaged Navier-Stokes equations, with the Boussinesq eddy-viscosity 
formulation, can be written in divergence form and then transformed from the Cartesian coordinate system 
(x, y, z) to the generalised system (£, The Jacobians of the coordinate transformation are defined as 
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where D is the determinant of the matrix J . Contravariant velocity components (£/, V, W) T = J 1 (u, v, w) T 
can be derived from the Cartesian components (u, u,tu) and then used in the transformed equations: 
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and p, P, and T are the fluid density, pressure, and temperature; and E is the total energy per unit volume. 
Pri = 0.72 is the laminar Prandtl number, Pr t = 0.90 is the turbulent Prandtl number, m is the laminar 
viscosity found from Sutherland’s law, and fit is the turbulent viscosity obtained from the k — e turbulence 
model. The total energy and pressure of a calorically perfect gas are related through the equation of state 
P = (7 - l)(P - pv • v/2), where 7 = 1.4 is the ratio of specific heat capacities. 

The eddy viscosity, fi ti is modelled as: 


= (2) 
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where jfc, the turbulence kinetic energy, and s, the dissipation rate of the turbulence kinetic energy, are 
obtained from the Launder and Spalding 2 high Reynolds number k — c turbulence model. C ^ is a scalar 
constant for isotropic turbulence and a wall function is used near solid boundaries. The k — e model has 
a number of benefits over the more widely used algebraic models. An explicit wake model is not required. 
This means that wakes do not have to be located and treated uniquely and are thus more likely to evolve and 
convect in a more realistic fashion. The wall function eliminates the need to sample information in a direction 
normal to solid boundaries thus making the calculation less sensitive to skewed or nonorthogonal mesh cells 
near solid walls. Furthermore this treatment, which makes no attempt to resolve the flow within the laminar 
viscous sublayer, requires fewer mesh cells and is thus less likely to be affected by the computational stiffness 
associated with highly stretched grids. Although it is true that heat transfer effects may not be modelled 
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accurately by a wall function, the stiffness of a low Reynolds number formulation and its corresponding poor 
convergence properties, are not present in the high Reynolds number model. The high Reynolds number 
k — e equations can be written: 
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for which the production rate of the turbulence kinetic energy can be defined as 
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The standard Launder and Spalding 1 constants 

= 0.09 h = 1.44 t 2 = 1.92 <x k = 1.0 <7, = 1.3 

were used throughout this work and no attempt was made to model nonisotropic turbulence. 

The transformed Navier-Stokes and k — e equations are discretized by a finite volume formulation that 
approximates the spatial differences as a net flux across the faces of each mesh cell. 3 Viscous stresses and 
geometric quantities are evaluated directly on the cell faces, while the flow and k — 5 variables are averaged 
over values found in adjacent cells. The unsteady equations are discretized into a linearized implicit delta 
form whose steady state solutions are not time-step dependent. To minimise the construction of the implicit 
operator, only the inviscid flux vectors are linearized. The flow equations can be written: 

[/ + t{S ( A + 6„B + 6'C)] A$,” fc = -At(S t (F - F u ) + S„(G - G v ) + S ( (H - #„)),"* (4) 


where A W n = W n + l - W n \ At is the time step size; 0 < & < 1 is a parameter governing the degree of 
implicitness; 6 and $ are cell- and face-centered central differences; / is the identity matrix; and A } B } and 
C are the inviscid flux Jacobian matrices relative to the vectors P t <5, and H . An implicit form of the k — € 
equations can also be written: 


+ pi&tfaAi" + S n B k9 + SfCkt - £*«)] AW k€ \ijk = ~A t(S^F kt + + $<H k4 - ( 5 ) 

where A kt1 B k€) and C kt are the inviscid flux Jacobian matrices relative to the inviscid terms found in the 
vectors F ktt G kti and R kt \ and E kt is the Jacobian matrix relative to the source vector § k4 . 
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Artificial Dissipation 

Explicit artificial dissipation terms must be added to this centered finite volume formulation to suppress 
possible odd and even point oscillations and shock overshoots. Fourth difference terms, similar to those used 
by Jameson 4 and Pulliam, 6 are added throughout the flow field to prevent odd-even decoupling, while second 
difference terms are added near shocks to stabilise the flow calculation. A local Mach number scaling is used 
to reduce the artificial dissipation in viscous-dominated flow regions and is similar to the treatment suggested 
by Flores and Holst The local Mach number scaling is normalised by the inflow Mach number and limited 
to a maximum value of unity. 

LU Approximate Factoriaation 

The equations are solved by approximately factoring the block-banded implicit operator into two block 
triangular ones. The LU factorisation, which is based on one-sided spatial differences, can be written for the 
flow' equations as: 

+ /i* At(6{ Ai + S~ B\ + £“Ci)j |V + A t(6^ A 2 + ^a)] = 

-A t(S<(? -P v ) + S„(d- G„)+S ( (i} - a v ) + f) n ijk (6) 

and for the k — e equations as follows: 
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where 8 + and 8~ are cell-centered forward and backward first differences; and A x = 0.5(A + fi\ A\I) and 
A 2 = 0.5(A - P\A\I) are the reconstructed flux Jacobian matrices, where |A| = max(|A^|) is the maximum 
absolute-valued eigenvalue of the Jacobian matrix A, fi « 1 is a scalar constant governing the amount of 
implicit dissipation produced by the matrix reconstructions, and I is the identity matrix. The vectors T 
and T kt represent the artificial dissipation terms added to the numerical schemes. These implicit systems 
of equations are solved by a two step procedure which, through back substitution, can be reduced to simple 
5x5 matrix systems at every mesh cell. These reduced systems can be written for the flow equations: 

1) Lower Sweep 
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2) Upper Sweep 
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and for the k — c equations 
1) Lower Sweep 
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2) Upper Sweep 
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(ii) 


Diagonally Inverted LU Factorisation 


The k — e equations are solved algebraically, while the flow equations are diagonally inverted using the 
similarity transformation: 1 
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and c is the local speed of sound. For a local time step defined as: 

A Cn 
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where Cn is the Courant number, the lower and upper sweeps, Eqs. 8 and 9 , can be approximated by the 
following scalar equations(m = 1, ..., 5 ). 
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2) Upper Sweep 
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where D ai = 0 , 5 (D a + P\A\I); D a , = 0 . 5 (Z> a — f)\A\I); and D a contains only the diagonal elements of 
the symmetric matrix Q~ l AQ [i.e., J ail = d a77 = d ajJ = U\ J a44 = U + c{$J 1 + £ y / 2 + £,^3); and 
5 = C c(£ z i 1 + £yJa ■+■ ^*^3)1 with similar terms for D Cl , and 

Initial and Boundary Conditions 

A uniform flow field based on an initial guess of the upstream conditions is used to start both the flow 
and k — e calculations. A no-slip condition is enforced along the solid boundaries and thus only pressures need 
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to be specified along the solid walls. These pressures are obtained from a three-dimensional interpretation 
of the normal momentum analysis developed by Risii. 7 

Total pressure, total temperature, and two flow angles are specified at subsonic inflow boundaries, while 
a one-dimensional Rieman invariant is extrapolated from the interior flow field. An inlet total temperature 
is fixed while a power law velocity profile is specified by varying the inlet total pressure. This treatment 
attempts to simulate an incoming boundary-layer profile. 8 The specified condition at a subsonic outflow 
boundary is a nonreflective treatment that attempts to minimise unwanted reflected waves. 0 Static pressures 
are obtained from a radial momentum analysis and then coupled with the incoming compatibility relation 
to produce a nonreflective boundary condition. 

An assumed two percent upstream turbulence intensity is used as an inflow boundary condition while 
the specified condition at the outflow boundary is a leroth order extrapolation of the k and c variables. 

Turbulent viscosities are set to iero along solid walls while a wall function 2 is used to evaluate the k — e 
variables in cells immediately adjacent to these solid boundaries. 

Steady State Calculations 

Local time-stepping and the multigrid method are used to accelerate the convergence of steady state 
calculations. 

A locally varying time step sixe, based on a constant Courant number, is defined throughout the flow 
field. This preconditioning creates a warped time integration that can accelerate calculations to a steady 
state without affecting the final solution. Local time stepping is used for both the flow and k — e equations, 
although the value of the Courant numbers used in their integration may vary. 

The flow algorithm is developed within the framework of the multigrid method to increase the efficiency 
of the time-marching procedure. Following the work of Jameson, 10 the basic flow algorithm is used to 
resolve the high frequency errors present on any current grid level (h), while the multigrid method is used to 
eliminate the lower frequencies through a sequence of calculations on coarser grids (2h i 4h i Sh } . . .). A four- 
level W-cycle, where a single flow calculation is performed on each grid level before transfering information 
to their respective coarser grids, was found to be the most efficient means of accelerating the steady state 
calculations. Coarse grid boundary conditions are identical to those on the fine grid, with the exception 
being that inflow/outflow conditions are updated only on the finest grid. The viscous fluxes are evaluated 
only on the finest grid, and thus influence the coarser grid calculations only through the transfered fine grid 
residual. 

The k — e equations are solved every multigrid cycle but only on the finest grid. They are not accelerated 
by the multigrid method defined above, which requires approximately 1.32 work units of computational effort. 
(One work unit is equivalent to a single Navier-Stokes calculation on the finest gnd.) There is little to be 
gained by multigridding the k-c equations between each fine grid flow solution since it is the time asymptote, 
rather than the solution at any current time level, that is ultimately required. Moreover, a coarse grid wall 
function, consistent with the fine grid flow field, is not easily defined and thus prohibits the simultaneous 
multigridding of both the flow and k - e equations. In addition, the viscous stresses cannot be adequately 
resolved on the coarser grids and thus their neglect at worst creates a high frequency error that, in any event, 
cannot be sustained by the multigrid method. 

Vector Sequence Acceleration Methods 

A number of vector sequence extrapolation methods have been used to accelerate steady state flow 
calculations. 1112,13 Most of these methods assume a slowly converging asymptotic linear behavior of the 
form 


W"* 1 = AW n +b (14) 

where W is the flow solution, 6 is an arbitrary vector, and A is the iteration matrix, whose largest eigenvalue 
has a magnitude close to unity. To remove the influence of the largest eigenvalue^), and thus increase 
the rate of convergence, a number of extrapolation methods have been developed. 1 One such method is 
Eddy’s 13 and Mesina’s 16 Reduced Rank Extrapolation. This method can be described as follows: 

Given the vector sequence W n \ where k < n, the following difference terms can be 

constructed 


R n = W n + l - W n 
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= R n+l - R n 

such that the solution of the normalised system of equations 

At— 1 

^ f (R n+i ,R n+3 )c, = -{R n+i ,R n ) i = o,i...,k-i (15) 

y=o 

(where (•/?"+ •', R*+>) is the inner product of the two vectors, R n+i and R n+1 ), can be used to construct an 
extrapolated solution: 


fc 

w ta = J2^ wn+1 ( 16 ) 
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where 
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7 3 = cy-1 ~ cy 1 < j < *“ 1 

Ik = Cfc-l 

For a fixed value of k t the accuracy of this extrapolated solution increases as n increases. In general it is 
true that the greatest amount of acceleration is produced when only a small value of k is needed to influence 
a small number of eigenvalues with large moduli. If a large number of dominant eigenvalues are of similar 
magnitude, then the extrapolation method may be of limited value since the manipulation of an equally 
large number of vectors would be required. As k increases, the solution of Eq. 15 becomes more sensitive to 
truncation errors and thus the accuracy of the extrapolated solution decreases. 

These vector sequence methods can also be used to analyse the convergence properties of any iterative 
method. 17 It can be shown 18 that, for n sufficiently large, the seros of the polynomial 

X>/A'=0 (17) 

y=o 

are estimates of the k largest eigenvalues of the iteration matrix A, Even though these values are only 
approximates, their locations with respect to one another can be used to analyse the convergence properties 
of the iterative method and predict the effectiveness of a vector sequence extrapolation. 

Results 

Numerical results for the flow through an annular turbine are used to illustrate this diagonally inverted 
LU implicit multigrid scheme’s ability to calculate three-dimensional compressible viscous flows. 

Turbomachinery calculations were performed on H-type grids consisting of 96x24x24 mesh cells in the 
throughflow, blade- to-blade, and radial directions, respectively. All calculations were performed on the 
NASA-Lewis CRAY X-MP, where approximately 1.5 million words of memory and 40 minutes of CPU time 
were required for a calculation consisting of 211 work units. 

The annular turbine, designed and tested at NASA Lewis, 19,20 is shown in Figure 1. This blade row is 
made up of 36 core turbine stator vanes, 38.10 mm high, with an axial chord of 38,23 mm. The stator has a 
tip diameter of 508 mm and a 0.85 hub-to-tip radius ratio. Mesh cells found immediately adjacent to solid 
walls are centered at distances 0.002 of an axial chord away, which correspond to a value of 60 < Y + < 200 
for the following flow calculations. 

Experimental test conditions of ambient axial inflow and a 0.65 hub-static to inlet-total pressure ratio 
produce a flow field with mean radius inlet and exit critical velocity ratios of 0.231 and 0.778, respectively. 
To match the upstream flow conditions (an inflow Mach number of 0.211), these nonrotating calculations 
were run with a 0.665 hub-static to inlet-total pressure ratio. 

The resulting flow field is fully subsonic and is compared with experimental data at three spanwise 
positions. Figures 2, 3, and 4 compare the calculated blade surface static pressure distributions (normaliied 
by inlet total pressure) at 13.3, 50, and 86.7 percent span with the experimented data produced by Goldman 
and Seasholtz. 19 Due to the bluntness of the blade’s trailing edge, the flow in this region is not sufficiently 
resolved, but the overall pressure distributions agree well with the experimental data. The calculations are 
most significant near the uncovered portion of the blade’s suction surface, where an inadequate turbulence 
model can greatly affect the accuracy of the flow solution. 


7 



Figure 5 shows the critical velocity ratio distributions along the 9, 50, and 81.7 percent radial span planes 
at the 155.8 percent axial chord location. The agreement between the computational and experimental results 
is extremely good and increases as one moves away from the hub and towards the tip region. The differences 
near the hub may be due to secondary flow and the non-isotropic turbulence found within this region. In 
general, the flow calculation has captured the viscous wake effects and agrees well with the experimental data. 
These results are obtained with the k - e turbulence model, which unlike the more widely used algebraic 
models, does not require the implimentation of an explicit wake model. 

Figure 6 shows the flow angles at the previously defined grid locations. The computational results 
compare extremely well with the experimental data at mid-span but deviate from the experimental data as 
one moves towards the hub and tip regions. These results are due to the underprediction of the boundary 
layer thicknesses and are not unexpected since the k — e model assumes fully turbulent flow and a transition 
model has not been incorporated into the calculation. 

Overall the numerical calculations compare extremely well with the experimental data and demonstrate 
the scheme’s ability to calculate three-dimensional viscous flows. 

To investigate the convergence properties of this multigrid scheme, a number of flow calculations, both 
viscous and inviscid, were performed. In all cases a Courant number of six was used to solve the flow 
equations while a Courant number of four was used for the k — € equations. The inviscid solutions were 
calculated on a 64x16x16 grid, which was of a relative fineness, comparable to that used for the viscous 
calculations. This was done in an attempt to insure that each of the flow calculations was performed on 
computational grids relevant to their respective length scales. It would seem to make little sense to perform 
viscous calculations on inviscid grids and vice versa. Thus, in order for meaningful comparisons to be made, 
each calculation, both viscous and inviscid, was performed on grids whose sizes were relevant to the physical 
problems being addressed. 

After four hundred work units on a single grid, the viscous flow calculation’s convergence history (Figure 
7) shows that the residual of the continuity equation has been reduced only two and a half orders of magni- 
tude. An abrupt ‘flattening out’ of the convergence history occurs after 300 work units and thus any further 
reduction would require significantly more iterations. The corresponding eigenvalue plot for this calculation 
is shown in Figure 8, where twenty-two of the largest eigenvalues are plotted in the complex plane. All of the 
eigenvalues have relatively equal magnitude and are thus located on a circle centered at the origin. This type 
of structure is characteristic of optimized SOR schemes, 2 * which are analogous to the LU implicit method in 
that both techniques require a back substitution step. Because of the large number of dominant eigenvalues, 
a vector sequence extrapolation would, in all likelihood, prove to be ineffective. The presence of virtually all 
of the eigenvalues would have to be accounted for before any amount of acceleration could be produced. 

The multigrid method is another way to accelerate convergence. Figure 9 shows the convergence history 
of the viscous multigrid calculation where a residual drop of approximately four and a half orders of magnitude 
was produced within 400 work units. The convergence acceleration is produced by a reduction in the 
magnitude of the largest eigenvalue(s), which are most likely associated with a low or intermediate spatial 
frequency. It is also interesting to discover that the multigrid method has altered the basic eigenvalue 
structure. The structure of the multigrid calculation, Figure 10, is clearly less clustered than that of the 
single grid solution. These results suggest that a multigrid scheme is more likely to benefit from a vector 
sequence extrapolation, an observation that has been confirmed in practice, by Reddy and Jacocks. 22 

These trends are also observed in the inviscid calculations. The convergence history (Figure 11) and 
eigenvalue structure (Figure 12) of the inviscid, single grid calculation again suggest the unlikelihood of 
using a vector sequence extrapolation method to accelerate convergence. The multigrid method produces a 
significant amount of convergence acceleration (Figure 13) and again has a less clustered eigenvalue structure 
(Figure 14). These calculations are more likely to be accelerated by a vector sequence extrapolation. Since 
both the inviscid and viscous eigenvalue plots are similar in structure (at least for the largest values), it 
would seem that, either inherently or by design, the convergence properties of the iterative procedures are 
inviscid ly dominated (inviscid in the sense that the high frequency modes most often associated with viscous 
effects do not seem to be dominating the convergence rates). This speculation is not unreasonable since it is 
also true that the viscous fluxes were not included in the implicit operator of the LU scheme. The differences 
in the amount of convergence acceleration produced by the viscous and inviscid multigrid calculations may 
be attributed more to differences in grid size and stretching, and less to viscous effects. 

Eriksson and Rizzi 23 found that the most persistant eigenmodes in their inviscid Runga-Kutta calcula- 
tions were most often associated with high spatial frequencies. They speculated that these high frequency 
modes could not be removed by the multigrid method, which operates on low frequency errors. Others have 
suggested that there would be little opportunity to accelerate viscous calculations since they presumably 
contain even more high frequency modes then do the inviscid equations. However, the existence of a num- 
ber of successful multigrid calculations indicates that there are significant amounts of low or intermediate 
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frequencies governing the convergence rates of a number of single grid calculations. The convergence accel- 
eration and modified eigenvalue structure produced by the multigrid method are evidence that a number of 
low frequency modes are dominant in these single grid calculations (Figures 7,8,11, and 12) and that they 
can be removed by the multigrid method (Figures 9,10,13 and 14). This analysis could also explain why a 
significant amount of acceleration occurs in the viscous calculations since presumably these low frequency 
modes exist in both the Euler and Navier-Stokes equations. 

In any event, the large number of dominant eigenvalues would prohibit the effectiveness of a vector 
sequence extrapolation in either the single or multigrid calculations. The eigenvalue plots suggest that at 
least 25-30 eigenvalues would have to be accounted for before any significant amount of acceleration could 
be achieved. This prediction is confirmed in Figure 15, where only a slight amount of acceleration (12 
percent more than the base multigrid convergence rate shown in Figure 9) is produced by the storage and 
manipulation of 25 viscous multigrid flow solutions. Although these results may seem discouraging, vector 
sequence methods are highly problem dependent and may be of greater value in other flow calculations. In 
addition, the eigenvalue estimates are extremely useful in the analysis of any iterative method. 

Concluding Remarks 

The convergence properties of an implicit multigrid scheme have been investigated for the calculation of 
complex three-dimensional flows. A vector extrapolation method was used to calculate eigenvalue estimates 
of the iterative procedure to gain a further understanding of the process by which the multigrid method 
accelerates calculations to steady state. 

Results illustrate both the scheme’s ability to calculate turbulent viscous flows and the convergence 
acceleration produced by the multigrid method. 
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